import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
import warnings

# 用于分析DPP4数据的各项统计学指标
warnings.filterwarnings('ignore')

# 读取CSV文件
df = pd.read_csv('/home/sunjun/ssd/dpp4/DPP-IV_RAQ.csv')

# 计算新列：DPP4_Digested减去No_DPP4
df['Difference'] = df['DPP4_Digested'] - df['No_DPP4']

# 样本数量
n = len(df['Difference'])
print(f'样本数量: {n}')

# 对于大样本，Shapiro-Wilk检验可能过于敏感，使用多种方法检验正态性
print('\n正态性检验结果:')

# 1. D'Agostino and Pearson's 检验
try:
    stat_dagostino, p_dagostino = stats.normaltest(df['Difference'])
    print(f"D'Agostino and Pearson's 检验: 统计量={stat_dagostino:.4f}, p值={p_dagostino:.4f}")
    if p_dagostino > 0.05:
        print(f"  根据D'Agostino检验，数据可能符合正态分布 (p > 0.05)")
    else:
        print(f"  根据D'Agostino检验，数据可能不符合正态分布 (p <= 0.05)")
except Exception as e:
    print(f"D'Agostino检验出错: {e}")

# 2. Anderson-Darling 检验
try:
    result = stats.anderson(df['Difference'], dist='norm')
    print(f"Anderson-Darling 检验: 统计量={result.statistic:.4f}")
    print("  临界值:")
    for i in range(len(result.critical_values)):
        sl, cv = result.significance_level[i], result.critical_values[i]
        if result.statistic < cv:
            print(f"    {sl}%: {cv:.4f} - 数据可能符合正态分布 (统计量 < 临界值)")
        else:
            print(f"    {sl}%: {cv:.4f} - 数据可能不符合正态分布 (统计量 >= 临界值)")
except Exception as e:
    print(f"Anderson-Darling检验出错: {e}")

# 3. Kolmogorov-Smirnov 检验
try:
    # 标准化数据
    data_standardized = (df['Difference'] - df['Difference'].mean()) / df['Difference'].std()
    stat_ks, p_ks = stats.kstest(data_standardized, 'norm')
    print(f"Kolmogorov-Smirnov 检验: 统计量={stat_ks:.4f}, p值={p_ks:.4f}")
    if p_ks > 0.05:
        print(f"  根据KS检验，数据可能符合正态分布 (p > 0.05)")
    else:
        print(f"  根据KS检验，数据可能不符合正态分布 (p <= 0.05)")
except Exception as e:
    print(f"KS检验出错: {e}")

# 4. Shapiro-Wilk 检验 (仅作参考，大样本时不太可靠)
if n <= 5000:  # Shapiro-Wilk通常只适用于较小样本
    try:
        stat_shapiro, p_shapiro = stats.shapiro(df['Difference'])
        print(f"Shapiro-Wilk 检验: 统计量={stat_shapiro:.4f}, p值={p_shapiro:.4f}")
        if p_shapiro > 0.05:
            print(f"  根据Shapiro-Wilk检验，数据可能符合正态分布 (p > 0.05)")
        else:
            print(f"  根据Shapiro-Wilk检验，数据可能不符合正态分布 (p <= 0.05)")
    except Exception as e:
        print(f"Shapiro-Wilk检验出错: {e}")
else:
    print("Shapiro-Wilk检验未执行 (样本量过大)")

# 绘制直方图和QQ图来可视化
plt.figure(figsize=(15, 10))

# 直方图与核密度估计
plt.subplot(2, 2, 1)
sns_plot = df['Difference'].plot(kind='hist', bins=30, alpha=0.7, color='skyblue', 
                                 edgecolor='black', density=True)
x = np.linspace(df['Difference'].min(), df['Difference'].max(), 100)
plt.plot(x, stats.norm.pdf(x, df['Difference'].mean(), df['Difference'].std()), 'r-', lw=2)
plt.title('Histogram with Normal Distribution Fit')
plt.xlabel('Difference (DPP4_Digested - No_DPP4)')
plt.ylabel('Density')

# QQ图
plt.subplot(2, 2, 2)
stats.probplot(df['Difference'], plot=plt)
plt.title('Q-Q Plot')

# 尝试拟合多种分布
print("\n拟合不同分布模型:")

# 使用scipy尝试几种常见分布
plt.subplot(2, 2, 3)
data_plot = plt.hist(df['Difference'], bins=30, alpha=0.5, density=True, label='Data')
x = np.linspace(df['Difference'].min(), df['Difference'].max(), 1000)

# 尝试几种常见分布
distributions = ['norm', 'laplace', 'logistic', 't', 'cauchy', 'gamma', 'lognorm']
colors = ['red', 'green', 'blue', 'purple', 'orange', 'brown', 'pink']
results = []

for i, dist_name in enumerate(distributions):
    try:
        dist = getattr(stats, dist_name)
        if dist_name == 't':
            params = dist.fit(df['Difference'], fscale=1)
        elif dist_name in ['gamma', 'lognorm']:
            # 对于某些分布，需要确保数据为正
            if dist_name == 'lognorm' and min(df['Difference']) <= 0:
                # 对于对数正态分布，数据必须为正
                shifted_data = df['Difference'] - min(df['Difference']) + 0.01
                params = dist.fit(shifted_data)
                print(f"{dist_name} 分布拟合: 数据已平移以满足正值要求")
            elif dist_name == 'gamma' and min(df['Difference']) <= 0:
                # 对于伽马分布，数据必须为正
                shifted_data = df['Difference'] - min(df['Difference']) + 0.01
                params = dist.fit(shifted_data)
                print(f"{dist_name} 分布拟合: 数据已平移以满足正值要求")
            else:
                params = dist.fit(df['Difference'])
        else:
            params = dist.fit(df['Difference'])
        
        y = dist.pdf(x, *params)
        plt.plot(x, y, colors[i], lw=2, label=f'{dist_name}')
        
        # 计算拟合优度 (使用KS检验)
        if dist_name in ['gamma', 'lognorm'] and min(df['Difference']) <= 0:
            # 对于已平移的数据，使用平移后的数据进行检验
            ks_stat, ks_pval = stats.kstest(shifted_data, dist_name, params)
        else:
            ks_stat, ks_pval = stats.kstest(df['Difference'], dist_name, params)
        
        print(f"{dist_name} 分布拟合: KS统计量={ks_stat:.4f}, p值={ks_pval:.4f}")
        
        # 存储结果用于排序
        results.append((dist_name, ks_stat, ks_pval, params))
    except Exception as e:
        print(f"{dist_name} 分布拟合出错: {e}")

# 按KS统计量排序（越小越好）
results.sort(key=lambda x: x[1])
print("\n分布拟合结果 (按KS统计量排序，越小越好):")
for i, (dist_name, ks_stat, ks_pval, params) in enumerate(results[:5]):
    print(f"{i+1}. {dist_name}: KS统计量={ks_stat:.4f}, p值={ks_pval:.4f}")
    print(f"   参数: {params}")

plt.legend()
plt.title('Distribution Fits (scipy)')
plt.xlabel('Difference (DPP4_Digested - No_DPP4)')
plt.ylabel('Density')

# 绘制累积分布函数比较
plt.subplot(2, 2, 4)
plt.hist(df['Difference'], bins=30, alpha=0.5, density=True, cumulative=True, label='Data CDF')

# 绘制前3个最佳拟合分布的CDF
best_dists = [result[0] for result in results[:3]]
for i, dist_name in enumerate(best_dists):
    if i < 3:  # 只绘制前3个
        dist = getattr(stats, dist_name)
        params = results[i][3]  # 获取参数
        
        if dist_name in ['gamma', 'lognorm'] and min(df['Difference']) <= 0:
            # 对于已平移的数据，需要在CDF计算中考虑平移
            shifted_x = x - min(df['Difference']) + 0.01
            y = dist.cdf(shifted_x, *params)
        else:
            y = dist.cdf(x, *params)
            
        plt.plot(x, y, colors[i], lw=2, label=f'{dist_name} CDF')

plt.legend()
plt.title('Cumulative Distribution Function Comparison')
plt.xlabel('Difference (DPP4_Digested - No_DPP4)')
plt.ylabel('Cumulative Probability')

plt.tight_layout()
plt.savefig('/home/sunjun/ssd/dpp4/distribution_analysis.png')

# 输出基本统计量
print('\n差值列的基本统计量:')
print(f'均值: {df["Difference"].mean():.4f}')
print(f'中位数: {df["Difference"].median():.4f}')
print(f'标准差: {df["Difference"].std():.4f}')
print(f'最小值: {df["Difference"].min():.4f}')
print(f'最大值: {df["Difference"].max():.4f}')
print(f'第一四分位数: {df["Difference"].quantile(0.25):.4f}')
print(f'第三四分位数: {df["Difference"].quantile(0.75):.4f}')
print(f'偏度: {stats.skew(df["Difference"]):.4f}')
print(f'峰度: {stats.kurtosis(df["Difference"]):.4f}')

# 保存带有差值列的数据到新文件
df.to_csv('/home/sunjun/ssd/dpp4/DPP-IV_test_with_difference.csv', index=False)

print('\n分析完成，结果已保存到 DPP-IV_test_with_difference.csv，图表已保存到 distribution_analysis.png')